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ABSTRACT 

A  brief  review  is  presented  of  various  problems  which  are  con¬ 
fronted  in  the  development  of  an  unsteady  finite  difference  poten¬ 
tial  code.  This  review  is  conducted  mainly  in  the  context  of  what 
is  done  for  a  typical  small  disturbance  and  full  potential  method. 
The  issues  discussed  Include  choice  of  equations,  linearization  and 
conservation,  differencing  schemes,  and  algorithm  development.  A 
number  of  applications,  including  unsteady  three-dimensional  rotor 
calculations,  are  demonstrated. 
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ONE  OF  THE  MOST  IMPORTANT  REALIZATIONS  in  fluid  mechanics  research  was 
that  transonic  flow,  for  all  its  apparent  complexity,  is  largely 
describable  by  potential  theory.  This  fact  was  clouded  by  a  thicket  of 
problems  concerning  tunnel  turbulence,  wall  and  scaling  effects,  and  sepa¬ 
ration.  We  now  know  that  many  of  these  problems  are  magnified  by  the 
inherent  susceptibility  of  the  inviscid  transonic  flow  to  unsteadiness. 

Of  course,  basic  flow  researchers  were  undoubtedly  prejudiced  in  favor  of 
potential  theory  from  the  outset  because  the  inviscid,  irrotational 
approximation  is  a  tremendous  simplification.  Nevertheless,  real  progress 
did  not  occur  until  the  explosive  development  of  finite  difference  methods 
which  occurred  in  the  1970* s. 

Spurred  on  by  the  promise  of  more  cruise-efficient  transport  air¬ 
craft,  our  greatest  progress  has  been  in  steady  flow  prediction.  However, 
unsteady  transonic  flows  are  important  because  the  aircraft  structural 
response  can  induce  large  and  heretofore  unpredictable  shock  excursions. 

The  flow  on  helicopter  rotor  tips  is  an  even  more  interesting  example  of 
unsteady  transonic  flow.  In  this  case,  unsteadiness  is  not  only  forced  by 
structural  deformation  but  also  by  a  free-stream  flow  which  rapidly  varies 
in  speed  and  direction  and  contains  large  wake  disturbances  from  previous 
blades.  And  so,  unsteady  transonic  flow  remains  a  rich  mine  of  important 
problems  waiting  to  be  solved.  Potential  methods  will  probably  play  a 
dominant  role  in  this  process  —  not  only  because  of  validity  but  also 
because  it  promises  to  be  the  most  efficient  method. 

Efficiency  is  a  more  important  matter  for  unsteady  computations  than 
for  the  steady  case.  This  is  because  we  require  the  resolution  of  physical 
time  and  do  not  have  the  benefit  of  acceleration  methods  that  are  used  in 
steady  problems.  In  general,  however,  the  unsteady  and  steady  methods  are 
much  the  same.  At  present,  the  most  versatile  and  efficient  unsteady 
methods  are  two-  and  three-dimensional  small  disturbance  codes,  but  full 
potential  methods  are  currently  under  rapid  development.  The  following 
discussion  will  review  many  of  the  important  algorithm  and  code  develop¬ 
ment  issues  in  the  context  of  small  disturbance  and  full  potential  methods. 

FORMULATIONS  OF  THE  PROBLEM 

The  starting  point  for  the  various  potential  formulations  is  the 
mass  conservation  and  Bernoulli's  equation  (shown  here  for  an  inertial 
reference  frame) 


P. 


Pt  +  V  •  (p V<|>)  »  0 


1 
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(1) 

(2) 


These  equations  have  the  advantage  (when  combined)  that  only  one  flow  vari¬ 
able,  need  be  solved  for.  However,  these  equations  are  only  an  approxi¬ 
mation  to  the  exact  inviscid  equations  because  they  assert  that  mass, 
energy  and  entropy  are  conserved  throughout  the  flow  field.  The  component 
of  momentum  normal  to  any  shock  is  not  conserved  in  this  approximation 
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(1,2)*-  This  is  a  notable  difference  from  the  Rankine-Hugoniot  equations 
where  shock  drag  comes  about  through  the  change  in  entropy  across  the 
shock.  The  error  thus  induced  is  generally  not  excessive  until  shock 
strengths  are  attained  that  involve  separation  and  the  necessity  to  abandon 
the  inviscid  approximation.  Conservative  formulations  are  required  to 
closely  approximate  Rankine-Hugoniot  results.  Nevertheless,  nonconserva¬ 
tive  formulations  are  still  actively  employed,  due  to  the  happy  accident 
that  errors  of  conservation  frequently  have  an  effect  similar  to  boundary- 
layer  corrections  on  shock  location. 

In  attacking  a  given  problem  it  is  first  necessary  to  express  Eqs.  (1) 
and  (2)  in  the  relevant  body-fixed  coordinate  system  for  the  problem  to  be 
solved.  The  transformation  (3)  which  encompasses  both  wings  or  rotors  in 
edgewise  motion  is 


r 1  *  U  t  +  R(t)r 
00 

t1  »  t 


(3) 


where  r*  »  (x^y',*1)  and  r  *  (x,y,z)  are  the  inertial  and  body-fixed 
coordinates  (see  Fig.  1),  respectively,  and  R(t)  is  the  rotation  matrix 

jcos  ft  t  -sin  71  t  0 

R(t)  -  sin  71  t  cos  71  t  0 

L  0  0  1 

(Note  that  for  no  rotation,  ft  »  0,  Eq.  (3)  reduces  to  the  usual  Galilean 
transformation.)  Under  this  transformation,  Eqs.  (1)  and  (2)  become 


-2-  - 
Poo 
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(4) 

(5) 


where  a  *  Uw(i  cosftt  -  j  sinftt)  -  ft  *  r  is  the  undisturbed  free-stream 
velocity  seen  by  an  observer  in  body-fixed  coordinates  and 

V  -  a  +  (6) 

is  the  flow  velocity  seen  in  the  moving  coordinates.  A  common  misunder¬ 
standing  concerns  the  validity  of  potential  theory  wi^ere  rotary  motion  is 
involved.  In  fact,  the  motion  of  the  coordinates,  -a,  has  no  bearing  on 
whether  or  not  a  potential  exists.  However,  since  a  need  not  be  an 
lrrotational  function,  V  is  not  generally  expressible  as  the  gradient  of 
a  full  potential.  Rather,  $  defines  a  disturbance  about  1,  as  seen  in 
Eq.  (6). 
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Numbers  in  parentheses  designate  References  at  end  of  paper. 


Equations  (4)  and  (5)  are  here  written  for  a  generalized,  moving  coor¬ 
dinate  system  where 


giving 

where 
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are  here  the  contravariant  velocity  components  and  have  the  form  of  Eq.  (6), 
J  is  the  Jacobian  | 3(5, n* O /3(x* ,y 1 ,z 1 ) |  and  V  is  the  cartesian 
gradient  operator.  In  the  same  manner,  Bernoulli's  equation  becomes  in 
these  coordinates 


P  *  1  - 


(Y-1) 


[?<♦,  +  Vc  +  Vn  +  Vt>  +  <Vs  +  Vn  +  V5>: 


*  (sy*5  +  Vn  +  ‘yV*  +  (V{  *  Vn  +  Vs’" 


3)  ^ 


(5a) 


where  we  use  the  nondimensionalizations  p  +  p/p^,  x  -►  x/i,  y  +  y/i, 
i  z/t%  $  +  and  T  *  tc^/i,  and  the  tilda  is  suppressed.  Here,  l  is 

a  reference  length  such  as  the  airfoil  chord. 

SMALL  DISTURBANCE  EQUATIONS — The  classical  small  disturbance  deriva¬ 
tion  involves  substituting  Eq.  (5)  into  Eq.  (4)  to  eliminate  p.  The 
resulting  equation  is  nondimensionalized,  scaled,  and  higher  order  terms 
are  eliminated,  subject  to  the  limit  process  that  (l-M2)/62/3  -  0(1)  as 
M  1  and  5+0.  The  resulting  equation  takes  the  form, 

A4>  +  B$  .  *  F  +  $  +  +  D$  (7) 

Ytt  rxt  x  rzz  yy  rxy 

where 

4  -*•  <t/u  162/3 
c 

D  *  frequency  of  unsteady  motion  (the  rota¬ 
tion  rate  for  a  rotor) 

M  *  U  /c  ,  a  characteristic  Mach  number 

U  ■  characteristic  speed  (V^  for  a  wing, 

c  DR  for  a  rotor) 

AR  »  R/i,  aspect  ratio 
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k  *  ft$,/Uc  reduced  frequency  (for  a  rotor 
k  -  1/AR) 

6  *  blade  thickness  ratio 
A  -  M2k2/S2/3 
B  -  2M2kf/62/3 
C  -  1/AR2S2/3 
D  -  Bg 

f  ■  y  +  ucos  t  (for  rotor),  or  1  (for  a  wing) 

F  ’  -157T  \  -  + 

g  *  x  4*  sin  t  (rotor),  or  1  (wing) 

V  »  ftR/Vw,  rotor  advance  ratio 
t  *  at* 
x  *  x?/£ 

y  «  y'/R  (R  is  either  a  rotor  radius  or  the 
wing  span  of  a  fixed  wing) 

z 1 6 1  /  3 
2  -  — — 

l  *  chord 


Equation  (7)  is  not  unique  and  an  assortment  of  modifications  and 
additions  have  been  made  to  improve  its  ability  to  handle  oblique  shocks 
(see  Ref.  4).  Nevertheless,  they  all  have  the  same  general  form  and  have 
essentially  identical  unsteady  terms. 

The  pressure  coefficient  in  the  small  disturbance  approximation  is 
given  by 


P  -  P 

c  -  T - 7  -  -  2ffi2/3  (K  +  k*T> 

t  pu  T 


(8) 


The  body  surface  boundary  condition  is  the  familiar  slope  condition 
applied  at  a  mean  surface. 


*  -  f(kb,  +  b  )  (9) 

Z  t  X 

where  the  body  surface  is  described  by  z  *  6b(x,y,t). 

There  is  an  additional  boundary  condition  that  pressure  be  continuous 
across  the  stagnation  streamline  behind  and  on  the  trailing  edge  of  an 
airfoil  (Kutta  condition).  Using  Eq.  (8),  this  condition  is  expressed  as 


Tx  +  krt  -  0  where  r 


(10) 


This  equation  expresses  the  downstream  convection  of  vorticity  (or  the 
potential  discontinuity).  It  happens  that  <j>  is  also  discontinuous  in 
the  airfoil  wake.  This  is  seen  by  applying  Eq?  (7)  (in  linearized  form) 
across  the  wake  and  substituting  Eq.  (10)  to  obtain 


where 


X-M'-[a+H8-t)]rxx-Cryy 
«  ■  ljM 
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This  quantity,  x*  is  zero  only  in  steady  two-dimensional  flow  or  steady 
nonlifting  three-dimensional  flow. 

Finally,  finite  difference  methods  require  far-field  boundary  condi¬ 
tions.  Typically,  one  specifies  some  combination  of  <p  »  0,  <p  ■  0,  or 
Cp  *  0  on  the  outer  boundaries  and  subsequently  relies  upon  anlarge 
boundary  distance  (often  over  1001)  to  dissipate  the  ensuing  wave  reflec¬ 
tions.  However,  it  has  been  shown  that  good  approximate  nonreflecting 
boundary  conditions  can  be  constructed  (5-8) .  For  example,  consider 
Eq.  (7)  in  its  two-dimensional,  linearized  form 


A*tt  +  B*xt  •  a<(,xx  + 


zz 


The  wave  information  for  a  plane  wave,  $  *  ei(u>t+£x+£z) ^  must  satisfy 

Aw2  +  *  C£2  +  Z2 

A  condition  which  prevents  reflection  of  upstream-moving  waves  at  the 
front  grid  boundary  is 


-B£  +  /(B^4AAaK^  -  4 An 2 


This  expression  does  not  transform  back  to  a  simple  local  differential 
operator.  However,  in  the  limit  n/£  +  0  (for  wavefronts  parallel  to  the 
upstream  boundary)  we  obtain  the  expression 

(  /B 2+4 Act  -  B\, 

V - 2A - )* 

whose  inverse  Fourier  transformation  yields 


/B2+4Aa  -  B 
- 2A 


(12) 


Similar  expressions  can  be  obtained  for  all  boundary  faces. 
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Although  its  limitations  are  many  and  well  known,  the  small  distur¬ 
bance  approximation  (with  judicious  use)  often  gives  excellent  results  and 
remains  in  active  use.  Perhaps  the  greatest  importance  of  small  distur¬ 
bance  theory  is  that  it  contains  (in  simplified  form)  all  the  important 
issues  contained  in  the  full  equations  and  hence  constitutes  a  good  test¬ 
ing  ground  for  future  techniques.  There  is,  perhaps,  one  exception  to  this 
statement.  It  concerns  the  fact  that  when  Eqs.  (4)  and  (5)  are  combined  to 
produce  a  single  equation  for  $,  a  nonconservative  formulation  results. 
This  problem  does  not  occur  in  classical  small  disturbance  theory. 

FULL  POTENTIAL  EQUATIONSH?or  computational  efficiency,  we  usually 
seek  to  formulate  the  potential  equations  so  that  we  need  only  one  depen¬ 
dent  variable.  This  is  easily  done,  since  Eq.  (4)  can  be  rewritten  as 


it  + 

3$  3t 


where 


1 32$  .  32*  ,  32*\ 
P\3x^"  +  3^2  +  dlz) 


<t  l£.l±+Q  &11+4  2e.ll 

x  3x  y  3$  By  z  3$  3z 


l£.m 

3* 


’2'Y(£+*x£  +  %^r+*z£) 


(13) 


is  a  none ommu ting  differential  operator.  For  this  equation  we  assume  the 
nonlinearizations  pertaining  to  Eq.  (5)  (for  the  previous  and  following 
equations  we  assume  pure  translational  motion) .  On  applying  the  isen- 
troplc  relations  we  obtain  the  familiar  textbook  form  of  the  potential 
equation 


*  +  2u<&  .  +  2v$  .  *  (c2-u2)$  +  (c2-v2)$  +  (c2-w2)<t  -  2$  $  $ 

tt  xt  yt  v  XX  yy  zz  x  y  xy 


—  2#  $  $  —  2$  $  $ 
x  z  xz  y  z  yz 


(14) 


The  problem  arises,  however,  that  equation  (14)  cannot  be  brought  back 
into  conservation  law  form  with  $  retained  as  the  dependent  variable. 

In  examining  Eq.  (4)  we  see  that  p  presents  no  problem  in  the 
special  derivative  terms  because  it  can  be  evaluated  from  the  previous 
time  step  (that  is,  it  is  taken  as  the  leading  term  in  a  Taylor  series  in 
time) .  The  required  implicit  special  terms  involving  $  then  come  about 
through  What  is  required  is  a  way  to  obtain  a  conservative  temporal 

function  of  $  from  the  term.  This  can  be  accomplished  through  the 

following  linearization, 

p  *  p0  +  (!t)0  <*-v  (15) 

where  the  subscript,  o,  denotes  a  nearby  known  state  or  solution.  With 
this  expansion  the  density  time  derivative  becomes 


which  provides  a  conservative  time-differenced  function  of  $  (9,10). 


(16) 
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The  above  linearization  can  be  avoided  for  a  conservative  formulation 
only  if  one  is  willing  to  retain  two  dependent  variables.  This  approach 
ban  lirfii  by  Ohlpm.iu  <m<l  Jameson  (11)  who  solve  tin*  llrnt-orilor 

ny*U  cm  t-oiiM t Ml  | tig  <»(  lln*  in.i •  im  ronriurvnt  Inn  i  Lou  anti  IWthou J  I  l  *  h  c,qii/.i- 
tion  rewritten  in  the  form 

*  0 

SPACIAL  DIFFERENCING 

The  primary  stability  issue  in  transonic  finite  difference  problems 
(steady  or  unsteady)  concerns  proper  spacial  differencing  in  the  subsonic 
and  supersonic  flow  regions.  Consider  first  the  simple  steady  two- 
dimensional  model  problem 

t1-*.2)  *XX  +  *yy  “  0  (17> 

and  let  the  difference  operators  VX,AX,V  ,  and  A  be  defined  as 

Vx<$>  =■  +  $£-i)/Ax,  Ax<$>  *  ($i+l  ~  $i)/Ax,  etc.  it  is  well  known  that  the 

following  difference  schemes 

(l-M2)  7A<j>  +  7A<(>  =  0  ,  M  <1  (18a) 

®  x  x  y  yT  « 


(l-M2)  V7^  +  7A^»0  ,  M  >  1  (18b) 

®  x  xT  y  yT  9  ®  v  7 

are,  respectively,  suitable  for  the  above  model  elliptic  and  hyperbolic 
problems.  The  differencing  equation  (Eq.  18b)  is  divergent  if  Moo  <  1 
while  (Eq.  18a)  is  convergent  for  >  1  only  if  ( Ax/[ (1-M»  )  Ay ] |  1, 

an  impractical  restriction  as  M*  1.  Murman  and  Cole  (12),  aware  of 
these  constraints,  introduced  type-dependent  differencing  for  the  transonic 
small  disturbance  equation.  In  this  approach  the  streamwlse  spatial  dif¬ 
ference  operator  is  either  central  or  backward,  depending  on  the  local 
Mach  number. 

In  order  to  obtain  a  stable  conservative  streamwise  difference  opera¬ 
tor  for  Eq.  (7),  it  is  first  necessary  to  time-linearize  the  nonlinear 
flux  term  F.  One  way  to  do  so  (13)  is 

This  time- linearized  flux  is  differenced  and  defined  at  each  midcell  of 
the  computational  grid  as 
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wnere  /  x 

Fl+is  “  -4x4n  ”  m2  \  2~)fAx'*,n  +  k(Y"1)7tl<,nJAxlt'n 

(jf)  -  a  -  M2  |(Y+l)fAx*n  +  k(Y-mt$n] 

and 

For  this  flux  operator  Murman’s  conservative  switching  scheme  (14)  is 

F*  *  V  1  ST  [£l('WFl->,)  +  (1-el-l)<Fi-^-Fi-3/2)]  <19’ 

where 

c.  -  f1  •  "l>0 

(o  ,  «4  <  0 

and 

Ma  -  a  -  M2(y+1)*x  -  M2(Y-l)k*t 

Note  that  for  greater  time  accuracy,  the  Crank-Nicolson  time  averaging  of 
spacial  operators^i  often  applied.  For  instance,  we  could  define  the 
flux  term  F  -  (F  +  Fn)/2  and  subsequently  apply  the  above  lineariza¬ 
tion  and  switching  scheme. 

The  nonconservative  full  potential  equation  (Eq.  14)  can  be  type- 
differenced  in  a  manner  very  similar  to  that  of  the  small  disturbance 
equation.  Equation  (14)  can  be  rearranged  in  the  cannonical  form  (15,16) 


♦tt  +  2q*st  *  <c2-«2)*ss  +  c2*nn  (20) 

where  s  is  a  coordinate  locally  aligned  to  the  stream  direction  and  n 
is  a  normal  coordinate.  For  the  two-dimensional  case 


j>  -  ~  (v2<j>  +2uv$  +u24>  ) 

nn  q2  xx  rxy  yy 


ss 


«  (u2<J>  +2  <J>  +v2<j>  ) 

q2  xx  uvTxy  Yyy 


and 


2q*st  *  2u*xt  +  2v*yt 


It  is  seen  that  the  steady  part  of  Eq.  (20)  has  the  same  form  as  Eq.  (18). 
Therefore  a  stable  scheme  results  when  $nn  is  always  central-differenced 
and  4>gg  is  type-differenced  based  on  the  sign  of  c2  -  q2. 

This  assertion  of  stability,  by  analogy  to  known  stable  methods,  can 
be  applied  to  the  full  potential  equation.  For  this  case  we  require  the 
following  difference  scheme  for  the  model  equation,  Eq.  (17). 
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VAd>-M2VV<t>  +  VA<fr»0  (21) 

X  X  oft  XX  y  y  \  / 

This  scheme,  though  less  accurate  than  Eq.  (18a),  is  stable  for  the 
entire  Mach  number  range  with  no  need  to  switch  difference  operators. 

Now  consider  the  spatial  derivative  term  3xp<f>x  of  Eq.  (4).  A 
local  linearization  of  this  term  gives 

3xp*x  "  3x [k)0  +  (♦xlf  +  p  h)o  (*“*o)]  (22) 

which  upon  substitution  of  Eq.  (13)  and  assumption  of  steady  small  distur¬ 
bance  flow  yields  the  approximate  expression, 

3x(pV  i  3x(poV  '  3x[(^l)0*x]  (23) 


If  we  further  assume  p  and  (u/c)  to  be  spacially  constant  we  have 

oo 

V“V>  ‘  “["x,-"2*,,,]  <« 

which  demonstrates  an  approximate  equivalence  between  the  full  potential 
equation  and  the  model  Eq.  (17).  The  importance  of  Eq.  (24)  is  that  the 
origin  of  the  second  on  the  right-hand  side  is  in  the  evaluation  of 

p  from  Bernoulli’s  equation.  Therefore,  the  differencing  scheme  of 
Eq.  (21)  is  significant  because  it  is  roughly  equivalent  to  evaluating 
3XP$X  with  a  centered  scheme  employing  an  upstream  biased  p.  The  sta¬ 
bility  of  this  density  biasing  has  been  demonstrated  in  Refs.  (17-19). 
This  density  biasing  will  be  expanded  upon  in  the  following  section. 

NUMERICAL  ALGORITHMS 


The  final  step  in  a  finite  difference  solution  scheme  is  the  effi¬ 
cient  solution  of  the  system  of  algebraic  equations  which  result  from 
the  various  discretizations.  Because  these  systems  are  all  too  large  to 
be  solved  efficiently  (in  spite  of  their  sparseness) ,  it  is  necessary  to 
reduce  them  to  a  more  manageable  form.  All  the  unsteady  schemes  today 
use  some  sort  of  approximate  factorization.  That  is,  the  system  matrix 
is  replaced  by  a  product  of  easily  solved  submatrices.  In  general,  this 
product  is  not  equal  to  the  original  system  matrix.  However,  it  is  pos¬ 
sible  to  keep  the  discrepancies  within  the  bounds  of  the  discretization 
error, 

SMALL  DISTURBANCE  EQUATION— The  germ  of  the  approximate  factorization 
idea  is  quite  old,  having  been  first  expressed  in  the  ADI  method.  A 
typical  ADI  scheme  for  Eq.  (7)  is 


Step  1.  ~  ?x(£-*n) 


D  F(i)  +  CV  A  An  +  V  A  d>n  +  DV 

v  r  \j  \r T  ZZ~ 


y  y 


vy*n 

Ay*0 


D  <  0 
D  >  0 
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Step  2.  ^  Vx($-$n) 


D  F($)  +  C7  4  +  V  i 

x  y  y  z  z 


n 


+  DV 

x 


Vy| 

Ayl 


D  <  0 
D  >  0 


Step  3.  A?tVt$n+1  +  B7t,?x<>n+1 


D  F(ijS)  +  CV  A  }  +  V  A 

X  y  y  2  Z 


♦n+1 


+  DV  < 
x 


D  <  0 
D  >  0 


Note  that  the  cross  derivative  term  is  treated  implicitly  and  is  always 
an  upwind  difference  (20,21),  The  advantage  of  this  ADI  scheme  is  that 
each  step  involves  a  simple  matrix  inversion  (steps  2  and  3  require  a 
tridiagonal  and  step  1  requires  a  quadradiagonal  inversion) .  Note  also 
that  these  steps  are  written  so  as  to  maximize  their  similarity  to  Eq.  (7). 
However,  they  are  never  actually  solved  in  this  form.  It  is  more  effici¬ 
ent  to  subtract  steps  1  and  2  from  steps  2  and  3,  respectively,  and  solve 
the  resulting  equations.  The  insertion  of  the  discretization  in 

the  above  algorithm  is  very  convenient  and  was  first  introduced  by 
Rizetta  and  Chen  (22). 

It  can  be  shown  that  the  above  ADI  scheme  constitutes  a  product  of 
terms  which  are  a  good  approximation  to  the  original  difference  equation. 
This  scheme  was  originally  devised  with  the  idea  that  each  step  should  be 
a  consistent  approximation  of  the  original  difference  equation.  However, 
if  this  constraint  is  relaxed  there  are  a  number  of  factorizations  which 
can  be  found.  To  see  how  this  is  done,  consider  the  model  equation 


which  is  differenced  as 


t  -  Bd  +  6 
*xt  M*xx  yyy 


7tVn+1  "  +  VyAy^n+l  <f°r  6  >  0) 


and  then  rearranged  as 


(7  -67  A  -7  A  )<j>n+1  -  7  4° 
x  x  x  y  y  x 


The  left-hand  side  of  this  equation  can  be  approximated  by  the  product 
(1-BAx) (Vx-VyAy)  if  the  error  term  BVxAyVy  is  eliminated.  This  elimi¬ 
nation  is  easily  accomplished  using  the  known  information  at  step  n. 

And  so,  an  appropriate  approximate  factorization  (commonly  called  the  AF2 
scheme)  is 


(1-67  )(7  -7  A  )<|>n+1  -  (7  +6A  7  A  )4>n 
x  x  yy  '  x  x  y  y  T 

This  gives  rise  to  a  two-step  procedure  —  one  step  being  a  bi-diagonal 
and  the  other  a  tri-diagonal  inversion.  A  further  simplification  is  to 
solve  for  ($n+i  -  $n)  rather  than  in  order  to  simplify  the  evalua¬ 

tion  of  the  right-hand  side. 
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CONSERVATIVE  FULL  POTENTIAL  METHOD-The  above  factorization  approach 
will  now  be  applied  to  the  unsteady  full  potential  equation.  We  will 
first  derive  a  simplified  form  of  the  equation  for  two  dimensions.  Con¬ 
sider  Eqs.  (4)  and  (5)  in  their  two-dimensional  form  subject  to  the 
mapping  5  ■  5(x),n  *  n(x)  and  assume  a  steady  free-stream  motion  allow¬ 
ing  $  to  be  expressed  as  $  >  M^x  +  $.  This  gives 

V»«> +  3e("T"  h)  *  'It-  0 ' 0  (25a) 

1 

P  -  jl  -  (y-1)[*t  +  \  (Cx2  «2  +  ny2  *y2  -M,2  )]j  (25b) 

On  application  of  the  linearization  of  Eq.  (15),  Eq.  (25a)  is  differenced 
as 


Vtpn(7t+5x2*e65  +  nz2^n<5n)(<>n+1”<>n)]  " 


+  6  (8°n  2  <S  *n+1)  +  v  Pn 

n  z  n  t 


(26) 


where  6  =  p2-'*’,  *  implies  division  by  J,  and  the  reference  state  0  is 
now  taken  to  be  the  previous  time  step  n.  The  spacial  operators  in 
Eq.  (26)  are  defined  as 


c  -  *1+1  ~  (  *1-1 

h  ~  2 

(  *i+l  “  (  *1-1 


6 «♦> 


(i+e)P.  +  (i-e)p .  .  “I 
*»ltl  - 4 - J - — J<W>1> 


(i+e)p, 


(27a) 

(27b) 


(27c) 
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5  (pn  2  5$) 
n  z  n 


j±i  1  PJ 


2  (Vi~V 


P1  +  pi-i 
*  ■  (a  —  <fc 

2  9j-l 


) 


(27d) 


Here,  AC  *  An  *  1  and  only  the  varying  indices  are  indicated.  The  param¬ 
eter  9  -  1  or  2  for  first-  or  second-order  spatial  accuracy  in  super¬ 
sonic  regions.  The  switching  parameter  v  is  defined  in  a  way  similar  to 
(17)  as 

v  -  [l  -  (p/p*)2]c  1  <_  c  <  10  | 

v  =  0  if  v  <  0  (i.e.,  subsonic)  >  (28) 

v  =  1  if  v  >  1  (i.e.,  supersonic)  I 


where  p*  is  density  evaluated  at  sonic  conditions. 

The  parameter  v  can  be  set  to  1  throughout,  but  accuracy  will  be 
impaired  unless  9  is  also  set  to  2.  The  operators  in  Eqs.  (27a)  and 
(27b)  assume  that  the  flow  will  be  supersonic  only  in  the  positive 
x-direction.  The  density  is  found  from  the  Bernoulli  equation  with 
(AC  «  An  -  1) : 


Ki+1  *i-l  *  , 

— - Vi 


Im 


Vj 


un+1  -  .♦“)  _ ,  ,n+i 
At  V 


The  metrics  C  and  n  are  obtained  from 
x  z 


Cx  “  xi+1  -  xi_1 


n„  - 


y  zj+i  "  zJ-i 


while  the  term  (£  2/J).,,  in  Eq.  (26)  is  formed  either  as 

X  l-r*5 


(5x2/J)i+1  +  Ux/J)i 


2 


(29a) 
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or 


VJ)i+i  *  (VJ)i 

2(xi+l-Xi) 


(29b) 


The  terms  /J)i-l/2  *  (nz2/J)  j+1/2  >  and  (nz2/J)j-i/2  receive  similar 
treatment.  If  Eq.  (19a)  is  used,  it  is  essential  to  add  -6^(p<i>53^ /J) 
to  the  right-hand  side  of  Eq.  (16)  to  subtract  out  a  numerical  truncation 
error  due  to  incomplete  metric  cancellation. 

Eq.  26  is  now  rearranged  into  delta_form;  it  is  now  to  be  solved  for 
*  4>n+l  .  $n#  For  example,  the  term  6^(5^  can  be  rearranged 

as 


VCx2  0D)VA*)  +  V?x2 

When  the  equation  is  put  into  delta  form  all  the  unknown  terms  are 
arranged  together,  resulting  in 

((0n/jn+1)[sx  +  Ux2)n+1  *"6^  +  <ny2)n+1  %n«n]-  h55(^x2/J)n+1 

-  h6n(ny2/J)n+1  pn6n)(<frn+1  -  *n) 

-  (6“*l/Ja)[sT  +  (Cx2)%“_\  +  (ny2)%rl5n>°-^l)  (30) 

+  (8n  -  an_1)  +  h(6^(c|/J)n+1  pn5c^R 
*  Sn<ny/J)'H'1  "■.V") 

which  can  be  approximately  factored  into  the  form 

{I  +  At(n  2)n+1  *"6  -  At(jn+1/0n)h6  (n  ^/J)1*1  Pn«.}x 

y  n  n  n  y  n 

(I  +  At(Cx2)“+1  -  At(Jll+1/8n)h«c(ex2/J)n+1  pn6c)(^n+1  -<^n) 

-  [1  +  (8n_1/en)(Jtt+1/Jn)]  (4>n-  ♦n“1)  -  (6n_1/6n)(Jn+1/Jn)(<j.n"1  -  4>n"2) 

(31) 

+  Ae<0n"l/8n)<J*fI/Jn)  [<€x2)“#J"l«5  +  (ny)D  -  4>n_1) 

+  At(Jn+1/Bn)((0n  -  pn_1)  +  h6  (52/J)n+1pn6_4.n 

**  X  s 

+  h60(r|2/j)"+1o“s^n) 
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This  equation  has  the  form 

L  L_(<J>n+1  -  $n)  ■  R  (28) 

n  S 

and  is  implemented  as  an  algorithm  as 


L  A<J>* 
n 

-  R 

(29) 

-  A<J>* 

n+1 

.  n  ,  .  ,  n 

♦ 

■  <t>  +  A<J> 

This  algorithm  requires  only  a  series  of  scalar  tridiagonal  inversions 
and  it  is  therefore  very  efficiently  implemented.  Computer  storage  equiv¬ 
alent  to  four  levels  of  $  have  to  be  supplied  with  p  computed  from  the 
Bernoulli  equation  as  needed. 

RESULTS  AND  APPLICATIONS 

Although  the  various  finite  difference  issues  are  very  similar  for 
the  unsteady  small  disturbance  and  full  potential  methods,  there  is  a 
large  difference  in  the  extent  to  which  these  methods  have  been  developed 
and  applied.  Naturally,  this  difference  reflects  the  relative  complexity 
of  the  two  methods. 

One  of  the  most  interesting  test  computations  demonstrating  the  use 
of  unsteady  potential  finite  difference  methods  concerns  the  flow  about 
a  pulsating  airfoil  (that  is,  one  having  a  time-varying  thickness). 
Although  the  idea  of  a  pulsating  airfoil  seems  farfetched,  it  can  be 
shown  by  considering  the  two-dimensional  small  disturbance  equation  to 
emulate  a  time-varying  free-stream  Mach  number.  And  free-stream  Mach 
number-variation  does  occur  for  an  advancing  helicopter  rotor.  Now  con¬ 
sider  a  parabolic  arc  airfoil  whose  mid-chord  thickness  varies  as 


T(t) 


f0.1[  10  -  t  +  6(t/15)2] (t/15)  3  ,  0  ^  t  <  15 

0.l[l0  -  15(22j*)  -  S^)2]^)3  ,  15  It  <  30 
Ko  ,  t  >  30 


where  t  -  t’UooM  and  is  nondimens ionalized  by  chord  lengths  traveled. 
Since  the  variation  takes  place  over  many  chords  of  travel  it  is  suitable 
to  invoke  the  low  frequency  approximation  for  the  small  disturbance  equa¬ 
tion  —  that  is,  all  time  derivatives  except  4>xt  are  ignored  in  Eq.  (7). 
(Physically,  ignoring  <f>tt  in  Eq.  (7)  amounts  to  assuming  an  infinite 
downstream  propagation  rate.)  Figures  2  and  3  show  a  comparison  of  the 
resulting  flow  computed  by  the  full  potential  method  (Goorjian,  Ref.  9) 
and  the  low  frequency  small  disturbance  equation  (Ballhaus  and  Steger, 
Ref.  23).  When  the  airfoil  is  thinning  it  becomes  subcritical  by  propa¬ 
gating  the  9hock  upstream  from  the  leading  edge.  The  two  approaches 


give  essentially  the  same  result.  The  most  significant  difference  between 
the  two  computations  is  perhaps  the  greater  dissipation  of  the  forward 
propagating  shock  for  the  full  potential  case.  This  can  be  controlled  by 
the  choice  of  upstream  density  biasing  function.  For  this  case,  density 
switching  was  employed.  Recall  that  it  is  possible  to  use  an  unswitched 
upstream  density  as  long  as  a  higher  order  difference  is  used  to  maintain 
accuracy.  This  is  demonstrated  in  Fig.  4  which  compares  switched  and 
unswitched  density  biasing  for  the  computation  of  the  steady  flow  on  a 
biconvex  profile. 

The  previously  mentioned  shock  motion  is  so  different  from  our 
previous  experience  that  it  demands  some  kind  of  experimental  study.  By 
coincidence,  such  a  study  was  performed  by  Tij deman  and  his  associates 
(24,25)  at  NLR,  the  Netherlands,  at  the  same  time  that  the  first  of  these 
computations  was  being  done.  They  acquired  detailed  flow  visualization 
and  loading  data  for  a  NACA64A006  airfoil  with  an  oscillating  flap.  They 
delineated  three  basic  types  of  shock  motion  caused  by  the  oscillating 
flap.  These  are: 

1.  Type  A.  The  shock  moves  nearly  sinusoidally  (only  the  lowest 
harmonic  was  measured)  with  a  phase  shift  relative  to  the  flap 
motion.  The  shock  strength  varies,  being  a  minimum  while  moving 
downstream  and  a  maximum  while  moving  upstream. 

2.  Type  B.  This  case  is  similar  to  Type  A  except  that  the  shock 
strength  variation  disappears  during  the  downstream  moving 
portion  of  its  cycle. 

3.  Type  C.  At  slightly  supercritical  conditions  there  is  no  super¬ 
sonic  region  for  a  large  portion  of  the  flap  cycle.  In  this 
case  the  airfoil  becomes  subcritical  by  propagating  the  shock 
upstream  off  of  the  leading  edge.  There  is  no  downstream  shock 
motion. 

The  ability  of  potential  methods  to  compute  these  shock  motions  was 
demonstrated  by  Ballhaus  and  Goorjian  (26).  In  this  work,  two-dimensional 
finite  difference  solutions  of  the  low  frequency  small  disturbance  equa¬ 
tion  were  obtained  by  the  ADI  approach  (program,  LTRAN2) .  The  computed 
Type  B  motion  is  shown  in  Fig.  5.  The  shock-motion  disappearance  and 
reappearance  are  shown  here  for  LTRAN2  and  the  Magnus-Yoshihara  code  (an 
isentropic  Euler  code).  Type  C  motion  is  illustrated  in  Fig.  6.  In  the 
succession  of  plots  shown,  the  shock  is  seen  to  form  at  mid-chord  and 
move  upstream  until  it  disappears  at  the  leading  edge.  The  shock  is  seen 
to  disappear  here  rather  than  visibly  propagate  upstream  as  in  Fig.  3, 
probably  due  to  numerical  dissipation.  It  is  interesting  that  the  condi¬ 
tions  for  which  these  various  types  of  shock  motion  were  computed  do  not 
compare  well  with  the  actual  experimental  conditions.  It  is  generally 
felt  now  (on  the  basis  of  comparing  different  codes  and  of  computations 
of  wall  effects)  that  the  discrepancy  is  due  to  tunnel-wall  and  possibly 
viscous  effects  rather  than  numerical  problems.  Nevertheless,  the  fact 
that  experiment  and  computation  produce  the  same  kinds  of  phenomena 
greatly  increases  the  significance  of  both. 

A  primary  application  of  unsteady  transonic  potential  methods  is  in 
the  prediction  of  loads  on  helicopter  rotors  in  forward  flight.  Although 
aeroelastic  effects  are  important,  in  this  case  the  main  source  of 
unsteadiness  is  in  the  flow  itself.  The  most  notable  distinction  between 
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the  fixed  and  rotary  wing  is  that  for  the  latter  the  free  stream  (in  body- 
fixed  coordinates)  is  constantly  accelerating  and  decelerating.  That  is, 
the  rotor  "sees"  the  free-stream  Mach  number  to  be  periodically  varying. 

The  reduced  frequency  of  this  variation  is  the  inverse  of  the  blade 
aspect  ratio.  Since  aspect  ratio  is  of  the  order  of  10  for  a  rotor,  it  is 
possible  to  invoke  the  low  frequency  approximation.  The  effectiveness  of 
the  small  disturbance  potential  methods  has  been  demonstrated  by  comparison 
of  computed  flows  with  measured  surface  pressure  on  a  nonlifting  rotor 
blade.  Figures  7  and  8  demonstrate  this  comparison  at  two  blade  azimuths, 

«  60°  and  120°  ,  respectively  (see  Fig.  7  for  definition  of  azimuth 
angle).  Measurements  and  computations  are  in  excellent  agreement.  The 
effects  of  unsteadiness  are  evident  in  this  figure  because  for  the  two 
azimuths  shown  the  chordwise  Mach  numbers  are  identical.  Yet  the  pres¬ 
sures  are  quite  different.  There  are  no  shocks  seen  at  *  60°  ,  but  they 
are  evident  at  120°  .  This  azimuthal  shock  asymmetry  is  not  explainable 
by  cross-flow  effects,  because  the  inboard  station  is  too  far  from  the 
tip.  Furthermore,  the  inboard  results  can  be  readily  obtained  by  two- 
dimensional  computations  (27) .  Very  often,  the  shock  motion  in  these 
rotor  computations  appears  to  be  Type  C,  the  upstream  propagating  type. 

This  is  seen  in  Fig.  9  (27)  which  shows  a  low  frequency  small  disturbance 
two-dimensional  computation  of  a  lifting  rotor  flow.  In  this  case  the 
blade  is  oscillating  and  also  sees  a  varying  free-stream  Mach  number. 

The  flow  on  the  bottom  surface  of  the  airfoil  is  seen  to  return  to  sub- 
critical  conditions  by  propagating  the  shock  upstream. 

The  most  fascinating  feature  of  the  above  computations  is  that  they 
are  so  easy  to  perform  compared  with  wind-tunnel  testing.  Yet  they  have 
often  proven  to  contain  most  of  the  essential  physics.  In  the  following 
discussion  we  shall  demonstrate  the  use  of  a  potential  computation  as  a 
"numerical  wind  tunnel"  to  explore  a  little-studied  but  possibly  very 
important  problem. 

An  unusual  unsteady  flow  feature  of  helicopter  rotors  is  that  they 
are  never  very  far  from  the  tip  vortex  of  a  preceding  blade  and  close 
blade/vortex  interactions  often  occur.  Under  certain  conditions  a 
blade  can  encounter  a  vortex  which  is  nearly  parallel  to  itself.  Such 
encounters  are  an  important  noise  source  and  require  modeling.  An 
extremely  simple  model  is  provided  by  a  two-dimensional  small  distur¬ 
bance  computation  of  a  near-vortex  encounter.  The  time  scale  for  this 
problem  (Z/U^)  is  very  brief  and  it  is  necessary  to  include  all  time 
derivatives.  The  vortex  is  introduced  as  the  edge  of  a  potential  dis¬ 
continuity  sheet  (Fig.  10)  which  is  stepped  through  the  computational 
grid.  The  vortex  is  moved  through  the  grid  in  a  prescribed  straight  line 
at  the  undisturbed  flow  speed.  The  strength  of  the  vortex  is  given  as 
an  effective  lift  coefficient,  Clv*  of  an  airfoil  having  the  same  circu¬ 
lation  as  the  vortex.  Note  that  while  the  sheet  describing  the  vortex 
in  Fig.  10  is  horizontal,  its  direction  is  irrelevant.  The  effect  of 
unsteadiness  in  these  computations  is  shown  in  Fig.  11,  which  compares 
blade  surface  pressure  distributions  for  a  fixed  and  moving  vortex  (the 
vortex,  located  at  mid-chord  and  0.96  chords  below  the  blade,  has  the 
strength,  Clv  *  0.1).  There  is  no  apparent  disturbance  for  the  unsteady 
moving  vortex  case  but  a  large  disturbance  for  the  steady  case.  This 
result  undoubtedly  reflects  the  fact  that  the  vortex  exerts  no  force  on 
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the  fluid  in  the  unsteady  case.  In  fact,  in  order  to  get  a  sizeable 
effect  on  the  surface  pressures  (other  than  the  expected  angle-of-attack 
variation)  it  is  necessary  to  bring  the  vortex  quite  close  to  the  blade 
and  Increase  its  strength.  Figure  12  shows  the  lift  variation  for  a 
blade  with  a  vortex  of  strength  Clv  *  0.4  whose  path  lies  0.25  chords 
beneath  the  blade.  As  the  vortex  passes  it  induces  a  shock  on  the  bottom 
surface.  This  shock  disappears  very  suddenly  with  no  obvious  upstream 
propagation.  Through  all  this  activity  the  upper  surface  seems  curiously 
unaffected. 

The  importance  of  such  a  computational  experiment  is  difficult  to 
assess  because  many  liberties  have  been  taken.  Nevertheless,  it  is  clear 
that  the  neglect  of  unsteadiness  would  be  to  greatly  overestimate  the 
effects  of  these  phenomena.  For  sufficiently  strong  or  close  vortices, 
large  and  rapid  bottom  surface  disturbances  can  be  computed  and  these 
must  have  acoustic  significance.  However,  this  could  be  mitigated  by 
allowing  the  vortex  to  move  freely  rather  than  follow  a  fixed  path. 
Undoubtedly,  an  assessment  of  these  sorts  of  computations  cannot  be  made 
without  some  experiments.  Thus  it  seems  that  the  effect  of  the  "numerical 
wind  tunnel"  can  be  to  guide  the  use  of  and  increase  the  need  for  the 
physical  wind  tunnel . 

CONCLUDING  REMARKS 

This  paper  has  made  no  attempt  to  treat  unsteady  potential  finite 
difference  methods  exhaustively  or  in  great  detail.  Rather,  a  typical 
small  disturbance  method,  and  a  full  potential  method  have  been  discussed 
in  parallel  in  the  context  of  the  issues  (choice  of  equations  and  boundary 
conditions,  linearizations,  discretizations,  and  algorithms)  which  arise 
in  all  code  development.  One  can  discuss  a  typical  method  because  the 
choices  taken  in  code  development  are  quite  few.  For  instance,  the  choice 
in  the  equations  centers  around  whether  one  wishes  to  have  one  or  two 
dependent  variables.  Linearizations  and  discretizations  vary  little. 

All  methods  employ  some  sort  of  upstream  biasing  and  nearly  the  same  time 
discretizations.  Algorithm  development  always  comes  down  to  some  sort  of 
approximate  factorization.  In  short,  although  there  is  much  development 
work  to  be  done,  the  general  area  of  potential  finite  difference  methods 
is  beginning  to  mature  quite  nicely.  We  shall  surely  see  the  development 
of  several  practical  three-dimensional  unsteady  full  potential  codes 
within  the  coming  year  or  two. 

There  remains  one  vital  issue  which  has  not  been  covered  in  this  dis¬ 
cussion:  the  subject  of  grid  generation.  On  this  point  finite  difference 

methods  are  often  more  of  an  art  than  a  science.  The  problem  is  not  as 
much  to  generate  a  mesh  as  to  know  the  effect  of  this  mesh  on  a  particular 
problem  and  solution  method.  This  is  especially  difficult  in  the  unsteady 
case  where  we  inav  have  the  presence  of  moving  flow  features  (shocks  and 
vortices^  which  require  resolution.  It  is  ironic  that  while  gridding  is 
the  most  fundamental  feature  of  finite  difference  methods,  the  topic  of 
grids  remains  c o  be  organized  into  an  organic  and  systematic  entity.  The 
requirements  lor  large-scale  unsteady  computations  must  certainly  change 
this  situation,  because  they  will  require  fast,  automatic,  and  reliable 
grid  schemes. 


18 


REFERENCES 

1.  J.  L.  Steger  and  B.  S.  Baldwin,  "Shock  Waves  and  Drag  in  the 
Numerical  Calculation  of  Isentropic  Transonic  Flow,"  NASA  TN  D-6997, 

Oct.  1972. 

2.  J.  vander  Vooren  and  J.  W.  Sloof,  "On  Inviscid  Isentropic  Flow 
Models  Used  for  Finite  Difference  Calculations  of  Two-Dimensional 
Transonic  Flows  with  Embedded  Shocks  About  Airfoils,"  NLR  MP  73024U. 

3.  M.  P.  Isom,  "Unsteady  Subsonic  and  Transonic  Potential  Flow  Over 
Helicopter  Rotor  Blades,"  NASA  CR-2463,  Oct.  1974 w 

4.  J,  vander  Vooren,  J.  W.‘  Sloof ,  G.  H.  Huizing,  and  A.  van  Essen, 
"Remarks  on  the  Suitability  of  Various  Transonic  Small  Perturbation 
Equations  to  Describe  Three-Dimensional  Transonic  Flow,"  presented  at 
Proc.  Symposium  Transonicum  II,  1975. 

5.  D.  Kwak,  f,Non-Reflecting  Far-Field  Boundary  Conditions  for 
Unsteady  Transonic  Flow  Computation,"  AIAA-80-1393,  presented  at  the  AIAA 
Fluid  and  Plasma  Dynamics  Conference,  Snowmass,  CO.,  July  14-16,  1980. 

6.  B.  Engquist  and  A.  Majda,  "Absorbing  Boundary  Conditions  for  the 
Numerical  Simulation  of  Waves,"  Math,  Comp.,  Vol.  31,  No.  139,  July  1977, 
pp.  629-651. 

7.  B.  Engquist  and  A.  Majda,  "Radiation  Boundary  Conditions  for 
Acoustic  and  Elastic  Wave  Calculations,"  Comm.  Pure  and  Applied  Math., 

Vol.  32,  May  1979,  pp.  313-357. 

8.  B.  Engquist  and  A.  Majda,  "Numerical  Radiation  Boundary  Condi¬ 
tions  for  Unsteady  Transonic  Flow,"  Journal  of  Comp.  Phys.,  Vol.  40, 

April  1981,  pp.  91-103. 

9.  P.  M,  Goorjian,  "Implicit  Computations  of  Unsteady  Transonic  Flow 
Governed  by  the  Full  Potential  Equation  in  Conservation  Form,"  AIAA-80- 
0150,  presented  at  the  AIAA  18th  Aerospace  Sciences  Meeting,  Pasadena,  CA., 
Jan.  14-16,  1980. 

10.  J.  L.  Steger  and  F.  X.  Caradonna,  "A  Conservative  Implicit  Finite 
Difference  Algorithm  for  the  Unsteady  Transonic  Full  Potential  Equation," 
AIAA-80-1368,  presented  at  the  AIAA  13th  Fluid  and  Plasma  Dynamics  Con¬ 
ference,  Snowmass,  CO.,  July  14-16,  1980. 

11.  R.  Chipman  and  A.  Jameson,  "Alternating  Direction  Implicit 
Algorithm  for  Unsteady  Potential  Flow,"  AIAA  Journal,  Vol.  20,  No.  1, 

Jan.  1982. 

12.  E.  M.  Murman  and  J.  D.  Cole,  "Calculation  of  Plane  Steady  Transonic 
Flows,"  AIAA  Journal,  Vol.  9,  No.  1,  1971,  pp.  114-121. 

13.  F.  X.  Caradonna  and  J.  L.  Steger,"  Implicit  Potential  Methods  for 
the  Solution  of  Transonic  Rotor  Flows,"  ARO  report  80-3,  Proceedings  of 
the  1980  Army  Numerical  Analysis  and  Computers  Conference. 

14.  E.  M.  Murman  and  J.  D.  Cole,  "Inviscid  Drag  at  Transonic  Speeds," 
presented  at  the  AIAA  7th  Fluid  and  Plasma  Dynamics  Conference,  Palo  Alto, 
California,  July  1973,  pp.  27-40. 

15.  K.  Isogai,  "Calculation  of  Unsteady  Transonic  Flow  Over  Oscil¬ 
lating  Airfoils  Using  the  Full  Potential  Equation,"  AIAA-77-448,  presented 
at  the  AIAA  Dynamics  Specialist  Conference,  San  Diego,  CA.,  Mar.  24-25, 
1977. 

16.  I-Chung  Chang,  "Unsteady  Transonic  Flow  Past  Airfoils  in  Rigid 
Body  Motion,"  DOE/ER/03077-170,  Courant  Institute  of  Mathematical 
Sciences,  Mar.  1981. 


19 


17.  T.  L.  Holst  and  W.  F.  Ballhaus,  "Fast  Conservative  Schemes  for 
the  Full  Potential  Equation  Applied  to  Transonic  Flows,"  AIAA  Journal, 

Vol.  17,  No.  2,  Feb.  1979,  pp.  145-152. 

18.  A.  Jameson,  "Transonic  Potential  Flow  Calculations  Using  Con¬ 
servative  Form,"  AIM  Second  Computational  Fluid  Dynamics  Conference 
Proceedings,  June  1975,  pp.  148-155. 

19.  M.  Hafez,  J.  South,  and  E.  Murman,  "Artifical  Compressibility 
Methods  for  Numerical  Solutions  of  the  Transonic  Full  Potential  Equation," 
AIM  Journal,  Vol.  17,  Aug.  1979,  pp.  838-844. 

20.  J.  J.  Chattot,  "Calculation  of  Three-Dimensional  Unsteady 
Transonic  Flows  Past  Helicopter  Blades,"  NASA  Technical  Paper  1721, 

Oct.  1980. 

21.  J.  J.  Philippe  and  J.  J.  Chattot,  "Experimental  and  Theoretical 
Studies  on  Helicopter  Blade  Tips  at  ONERA,"  Paper  No.  46,  presented  at 
the  Sixth  European  Rotorcraft  and  Powered  Lift  Aircraft  Forum,  Bristol, 
U.K.,  Sept.  1980,  pp.  16-19. 

22.  D.  P.  Rizzetta  and  W.  C.  Chin,  "Effect  of  Frequency  in  Unsteady 
Transonic  Flow,"  AIM  Journal,  Vol.  17,  No.  7,  July  1979,  pp.  779-781. 

23.  W.  F.  Ballhaus  and  J.  L.  Steger,  "Implicit  Approximate- 
Factorization  Schemes  for  the  Low-Frequency  Transonic  Equation,"  NASA 
TM  X-73,082,  Nov.  1975. 

24.  H.  Tijdeman,  "Investigations  of  the  Transonic  Flow  Around  Oscil¬ 
lating  Airfoils,"  Ph.D.  Thesis,  National  Aerospace  Laboratory 
(Netherlands),  TR  77090,  Oct.  1977. 

25.  H.  Tijdeman,  P.  Schippers,  and  A.  J.  Persoon,  "Unsteady  Airloads 
on  an  Oscillating  Supercritical  Airfoil,"  AGARD  CP-226,  Unsteady  Air¬ 
loads  in  Separated  and  Transonic  Flow,  Apr.  1977. 

26.  W.  F.  Ballhaus  and  P.  M.  Goorjian,  "Implicit  Finite-Difference 
Computations  of  Unsteady  Transonic  Flows  About  Airfoils,  Including  the 
Treatment  of  Irregular  Shock-Wave  Motions,"  AIM  Journal,  Vol.  15, 

No.  12,  Dec.  1977,  pp.  1728-1735. 

27.  F.  X.  Caradonna  and  J.  J.  Philippe,  "The  Flow  Over  a  Helicopter 
Blade  Tip  in  the  Transonic  Regime,"  Vertica,  Vol.  2,  1978,  pp.  43-60. 


1  ***■%►- 


22 


FULL  POTENTIAL  EQUATION 
-  0.0625 


Fig*  3  -  Pressure  variation  on  a  pulsating  airfoil — airfoil  thinning 
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Fig.  11  -  Steady  and  unsteady  flow  computations  for  the 
blade/vortex  interaction  problem 
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